While it is common knowledge that the two-body Kepler problem can be solved exactly as an orbit equation, it is perhaps less commonly known that there is also a complete solution as a function of time in terms of Kapteyn series. These are series of the form

$\sum _{k=0}^{\infty }{c}_k{J}_k\left(kz\right)$

where the index of the Bessel function appears as a factor in its argument. Such series with appropriate coefficients will converge as long as the quantity

$\left|\frac{zexp\left(\sqrt{1-z^2}\right)}{1+\sqrt{1-z^2}}\right|$

remains finite, which is true even for extremely eccentric orbits in the Kepler problem. This presentation will evaluate the explicit temporal expansions for dynamic variables of interest and their combinations. It expands on the cursory treatment given in §17.21 of Watson’s venerable Bessel Functions.

The equation of an orbit with semimajor axis a and eccentricity e, along with its parametrization in terms of the eccentric anomaly u, is

$r=\frac{a(1-e^2)}{1+ecos\phi }=a(1-ecosu)$

The relation between the two angles in the invariant plane can be rearranged to give

$\frac{1-e^2}{1+ecos\phi }=1-ecosu\Rightarrow \begin{array}{l}cos\phi =\frac{cosu-e}{1-ecosu}\\ sin\phi =\frac{\sqrt{1-e^2}sinu}{1-ecosu}\end{array}$

from which one immediately has

$x=rcos\phi =a(cosu-e)y=rsin\phi =a\sqrt{1-e^2}sinu$

Note that the x-coordinate is related linearly to the radial variable. Integrating the energy equation in terms of the eccentric anomaly produces Kepler’s equation,

$\omega t=u-esinu\omega =\sqrt{\frac{k}{m{a}^3}}$

where the frequency of the orbit is dependent on the reduced mass m and the gravitational interaction constant k as well as the semimajor axis. Differentiating both sides of Kepler’s equation

$d\left(\omega t\right)\equiv d\tau =(1-ecosu)du$

provides a differential relationship between the variables in the equation. The symbol τ as defined here will be used as shorthand to simplify integral expressions.

The Kapteyn series representing dynamic variables arise from Fourier sine and cosine expansions. Given a periodic function of the eccentric anomaly, a sine series is determined by

$f\left(u\right)=\sum _{k=1}^{\infty }{s}_{k}sink\omega t{s}_k=\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}f\left(u\right)sink\tau d\tau$

while a cosine series is determined by

$f\left(u\right)=\sum _{k=0}^{\infty }{c}_{k}cosk\omega t\begin{array}{l}{c}_0=\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}f\left(u\right)d\tau \\ c_k=\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}f\left(u\right)cosk\tau d\tau \end{array}$

The choice between the two expansions depends on the symmetry of the periodic function: an odd function will be expanded as a sine series and an even function as a cosine series.

Evaluation of the coefficients necessitates an integral representation of Bessel functions. While there are various ways to write such a representation, the one that is most useful in this context is

$J_n\left(z\right)=\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}cos(n\theta -zsin\theta )d\theta$

which holds for n integral. As an aside, the corresponding integral with a sine instead of a cosine is related to the Bessel function $Y_n\left(z\right)$ .

Differentiation with respect to the argument z under the integral provides an integral representation of the Bessel function derivative:

$J_n^{\prime }\left(z\right)=\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}sin\theta sin(n\theta -zsin\theta )d\theta$

Application of product identities for circular functions gives

$\begin{array}{l}\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}cosm\theta cos(n\theta -zsin\theta )d\theta =\frac{1}{2}\left[J_{n-m}\right(z)+J_{n+m}(z\left)\right]\\ \frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}sinm\theta sin(n\theta -zsin\theta )d\theta =\frac{1}{2}\left[J_{n-m}\right(z)-J_{n+m}(z\left)\right]\end{array}$

Two recursion relations for Bessel functions

$\begin{array}{l}{J}_{n-1}\left(z\right)+J_{n+1}\left(z\right)=\frac{2n}{z}{J}_n\left(z\right)\\ J_{n-1}\left(z\right)-J_{n+1}\left(z\right)=2J_n^{\prime }\left(z\right)\end{array}$

will be useful in what follows. The second can be seen immediately in the integral representations above, while the first implies

$\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}cos\theta cos(n\theta -zsin\theta )d\theta =\frac{n}{z}{J}_n\left(z\right)$

Keep in mind that a prime will always indicate differentiation with respect to the entire argument of the Bessel function. A second derivative can always be replaced as

$J_n^{″}\left(z\right)=-\frac{1}{z}{J}_n^{\prime }\left(z\right)+[\frac{n^2}{z^2}-1]J_n\left(z\right)$

using the defining differential equation for Bessel functions.

Inversion of Kepler’s equation requires an expansion for the sine of the eccentric anomaly. The twist to the integral involved and those that follow is that the integration is with respect to the mean anomaly τ, not the eccentric anomaly. This will add a factor to the integrand of the differential relationship between the variables in the equation.

Integrals can always be evaluated by expanding products of circular functions and applying recursion relations to the resulting Bessel functions, but often an intermediate integration by parts is more efficient. This is particularly true with a single differential factor in the integrand.

The coefficients for a sine expansion of the sine of the eccentric anomaly are

$\begin{array}{l}\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}sinusink\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}sinusin\left[k(u-esinu)\right](1-ecosu)du\\ =\frac{2}{\pi k}\underset{0}{\overset{\pi }{\int }}cosucos\left[k(u-esinu)\right]du=\frac{2}{ke}{J}_k\left(ke\right)\end{array}$

where the endpoint contributions to the integral are zero. The expansion for the sine of the eccentric anomaly is thus

$sinu=\frac{2}{e}\sum _{k=1}^{\infty }\frac{J_k\left(ke\right)}{k}sink\omega t$

and that for the eccentric anomaly itself is

$u=\omega t+2\sum _{k=1}^{\infty }\frac{J_k\left(ke\right)}{k}sink\omega t$

The corresponding cosine expansion for the cosine of the eccentric anomaly has a constant term

$\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}cosud\tau =\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}cosu(1-ecosu)du=-\frac{e}{2}$

and the remaining coefficients are

$\begin{array}{l}\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}cosucosk\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}cosucos\left[k(u-esinu)\right](1-ecosu)du\\ =\frac{2}{\pi k}\underset{0}{\overset{\pi }{\int }}sinusin\left[k(u-esinu)\right]du=\frac{2}{k}{J}_k^{\prime }\left(ke\right)\end{array}$

where endpoint contributions to the integral are again zero. The complete expansion of the cosine of the eccentric anomaly is

$cosu=-\frac{e}{2}+2\sum _{k=1}^{\infty }\frac{J_k^{\prime }\left(ke\right)}{k}cosk\omega t$

Having expansions of both the sine and cosine of the eccentric anomaly immediately gives expansions for the radial variable and its Cartesian components. The expansion of the radial variable is

$\frac{r}{a}=1+\frac{e^2}{2}-2e\sum _{k=1}^{\infty }\frac{J_k^{\prime }\left(ke\right)}{k}cosk\omega t$

while those of its two components are

$\frac{x}{a}=-\frac{3e}{2}+2\sum _{k=1}^{\infty }\frac{J_k^{\prime }\left(ke\right)}{k}cosk\omega t$

$\frac{y}{a}=\frac{2\sqrt{1-e^2}}{e}\sum _{k=1}^{\infty }\frac{J_k\left(ke\right)}{k}sink\omega t$

An expansion of the inverse of the radial variable is particularly easy because of cancellations of linear factors. The constant term in its cosine expansion is

$\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}\frac{a}{r}d\tau =\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}du=1$

while the remaining coefficients are

$\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}\frac{a}{r}cosk\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}cos\left[k\right(u-esinu\left)\right]du=2J_k\left(ke\right)$

so that the complete expansion is

$\frac{a}{r}=1+2\sum _{k=1}^{\infty }{J}_k\left(ke\right)cosk\omega t$

The same cancellations occur in expansions for the sine and cosine of the polar angle in the invariant plane. The constant term in the cosine expansion is

$\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}cos\phi d\tau =\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}(cosu-e)du=-e$

and the remaining coefficients are

$\begin{array}{l}\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}cos\phi cosk\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}(cosu-e)cos\left[k\right(u-esinu\left)\right]du\\ =\frac{2(1-e^2)}{e}{J}_k\left(ke\right)\end{array}$

so that the complete expansion is

$cos\phi =-e+\frac{2(1-e^2)}{e}\sum _{k=1}^{\infty }{J}_k\left(ke\right)cosk\omega t$

This expression can also be found using the relationship

$cos\phi =-\frac{1}{e}+\frac{1-e^2}{e}\frac{a}{r}$

between dynamic variables and the previous result. The coefficients for the corresponding sine expansion are

$\begin{array}{l}\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}sin\phi sink\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}\sqrt{1-e^2}sinusin\left[k\right(u-esinu\left)\right]du\\ =2\sqrt{1-e^2}{J}_k^{\prime }\left(ke\right)\end{array}$

with the result

$sin\phi =2\sqrt{1-e^2}\sum _{k=1}^{\infty }{J}_k^{\prime }\left(ke\right)sink\omega t$

An explicit expansion of the polar angle φ itself is possible, but it is much more complicated than either of the last two expressions and will not be included.

With the basic dynamic variables in place, one can begin to consider expansions for some of their combinations. First consider a cosine expansion of the square of the sine of the eccentric anomaly. The constant term is

$\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}{sin}^{2}ud\tau =\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}{sin}^{2}u(1-ecosu)du=\frac{1}{2}$

and the remaining coefficients are

$\begin{array}{l}\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}{sin}^{2}ucosk\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}{sin}^{2}ucos\left[k(u-esinu)\right](1-ecosu)du\\ =-\frac{2}{\pi k}\underset{0}{\overset{\pi }{\int }}sin2usin\left[k(u-esinu)\right]du\\ =\frac{1}{k}[J_{k+2}\left(ke\right)-J_{k-2}\left(ke\right)]\\ =\frac{1}{k}[\frac{2(k+1)}{ke}{J}_{k+1}\left(ke\right)-\frac{2(k-1)}{ke}{J}_{k-1}\left(ke\right)]\\ =\frac{4}{k^2{e}^2}{J}_k\left(ke\right)-\frac{4}{ke}{J}_k^{\prime }\left(ke\right)\end{array}$

where $J_0\left(ke\right)$ is added and subtracted on the third line. The complete expansion for this square is thus

${sin}^{2}u=\frac{1}{2}+\frac{4}{e^2}\sum _{k=1}^{\infty }\frac{J_k\left(ke\right)}{k^2}cosk\omega t-\frac{4}{e}\sum _{k=1}^{\infty }\frac{J_k^{\prime }\left(ke\right)}{k}cosk\omega t$

The expansion for the square of the cosine of the eccentric anomaly follows from the simple relationship between the two circular functions:

${cos}^{2}u=\frac{1}{2}-\frac{4}{e^2}\sum _{k=1}^{\infty }\frac{J_k\left(ke\right)}{k^2}cosk\omega t+\frac{4}{e}\sum _{k=1}^{\infty }\frac{J_k^{\prime }\left(ke\right)}{k}cosk\omega t$

Another combination that follows from the integal in the last evaluation is a sine series for the product of the sine of the eccentricity and the cosine of the angular variable in the invariant plane. The coefficients are

$\begin{array}{l}\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}sinucos\phi sink\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}sinu(cosu-e)sin\left[k\right(u-esinu\left)\right]du\\ =-2e{J}_k^{\prime }\left(ke\right)+\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}sin2usin\left[k\right(u-esinu\left)\right]du\\ =-\frac{2}{k{e}^2}{J}_k\left(ke\right)+\frac{2(1-e^2)}{e}{J}_k^{\prime }\left(ke\right)\end{array}$

with the result

$sinucos\phi =-\frac{2}{e^2}\sum _{k=1}^{\infty }\frac{J_k\left(ke\right)}{k}sink\omega t+\frac{2(1-e^2)}{e}\sum _{k=1}^{\infty }{J}_k^{\prime }\left(ke\right)sink\omega t$

The product of the sine of the eccentricity and the sine of the angular variable is a cosine series whose constant term is

$\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}sinusin\phi d\tau =\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}\sqrt{1-e^2}{sin}^{2}udu=\frac{1}{2}\sqrt{1-e^2}$

The remaining coefficients are

$\begin{array}{l}\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}sinusin\phi cosk\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}\sqrt{1-e^2}{sin}^{2}ucos\left[k\right(u-esinu\left)\right]du\\ =\sqrt{1-e^2}{J}_k\left(ke\right)-\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}cos2ucos\left[k\right(u-esinu\left)\right]du\\ =\sqrt{1-e^2}[J_k\left(ke\right)-\frac{1}{2}{J}_{k+2}\left(ke\right)-\frac{1}{2}{J}_{k-2}\left(ke\right)]\\ =\sqrt{1-e^2}[2J_k\left(ke\right)-\frac{(k+1)}{ke}{J}_{k+1}\left(ke\right)-\frac{(k-1)}{ke}{J}_{k-1}\left(ke\right)]\\ =\sqrt{1-e^2}[-\frac{2(1-e^2)}{e^2}{J}_k\left(ke\right)+\frac{2}{ke}{J}_k^{\prime }\left(ke\right)]\end{array}$

so that the complete result is

$sinusin\phi =\sqrt{1-e^2}[\frac{1}{2}-\frac{2(1-e^2)}{e^2}\sum _{k=1}^{\infty }{J}_k\left(ke\right)cosk\omega t+\frac{2}{e}\sum _{k=1}^{\infty }\frac{J_k^{\prime }\left(ke\right)}{k}cosk\omega t]$

Finally, since the equations of motion in Cartesian coordinates are

$\begin{array}{l}\frac{d^{2}x}{d{t}^2}=-\frac{k}{m}\frac{x}{r^3}=-{\omega }^2{a}^3\frac{cos\phi }{r^2}\\ \frac{d^{2}y}{d{t}^2}=-\frac{k}{m}\frac{y}{r^3}=-{\omega }^2{a}^3\frac{sin\phi }{r^2}\end{array}$

one can differentiate existing expansions for the immediate results

$\frac{a^2}{r^2}cos\phi =2\sum _{k=1}^{\infty }k{J}_k^{\prime }\left(ke\right)cosk\omega t$

$\frac{a^2}{r^2}sin\phi =\frac{2\sqrt{1-e^2}}{e}\sum _{k=1}^{\infty }k{J}_k\left(ke\right)sink\omega t$

The direct evaluation of these expansions is most easily accomplished with an intermediate integration by parts using the integrals

$\begin{array}{l}\int du\frac{cosu-e}{(1-ecosu{)}^2}=\frac{sinu}{1-ecosu}\\ \int du\frac{sinu}{(1-ecosu{)}^2}=-\frac{1}{e(1-ecosu)}\end{array}$

where endpoint contributions to integrals are zero in both cases. The coefficients for the cosine expansion are

$\begin{array}{l}\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}\frac{a^2}{r^2}cos\phi cosk\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}\frac{cosu-e}{(1-ecosu{)}^2}cos\left[k\right(u-esinu\left)\right]du\\ =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}ksinusin\left[k\right(u-esinu\left)\right]du=2k{J}_k^{\prime }\left(ke\right)\end{array}$

while the coefficients for the sine series are

$\begin{array}{l}\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}\frac{a^2}{r^2}sin\phi sink\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}\frac{\sqrt{1-e^2}sinu}{(1-ecosu{)}^2}sin\left[k\right(u-esinu\left)\right]du\\ =\frac{2}{\pi }\frac{\sqrt{1-e^2}}{e}\underset{0}{\overset{\pi }{\int }}kcos\left[k\right(u-esinu\left)\right]du=\frac{2k\sqrt{1-e^2}}{e}{J}_k\left(ke\right)\end{array}$

Both are identical as expected with the results from the equations of motion. A third way to get to the cosine expansion is via partial differentiation of the expansion of the inverse radial variable:

$\begin{array}{l}\frac{\partial }{\partial e}[\frac{a}{r}=\frac{1}{1-ecosu}=1+2\sum _{k=1}^{\infty }{J}_k\left(ke\right)cosk\omega t]\\ \frac{1}{(1-ecosu{)}^2}\frac{\partial [ecosu]}{\partial e}=2\sum _{k=1}^{\infty }k{J}_k^{\prime }\left(ke\right)cosk\omega t\\ \\ \frac{\partial }{\partial e}[ecosu]=\frac{\partial }{\partial e}[-\frac{e^2}{2}+2e\sum _{k=1}^{\infty }\frac{J_k^{\prime }\left(ke\right)}{k}cosk\omega t]\\ =-e+2\sum _{k=1}^{\infty }\frac{J_k^{\prime }\left(ke\right)}{k}cosk\omega t+2e\sum _{k=1}^{\infty }{J}_k^{\prime \prime }\left(ke\right)cosk\omega t\\ =-e+\frac{2(1-e^2)}{e}\sum _{k=1}^{\infty }{J}_k\left(ke\right)=cos\phi \end{array}$

where defining differential equation for Bessel functions has been used to replace the second derivative. This last somewhat surprising evaluation raises the question of what relationships hold among the various Kapteyn series appearing in this presentation and the extent to which such relationships can produce further expansions.

To facilitate comparisons of expressions, introduce the notation for the series

$\begin{array}{l}{S}_n=\sum _{k=1}^{\infty }{k}^n{J}_k\left(ke\right)sink\tau \\ C_n=\sum _{k=1}^{\infty }{k}^n{J}_k\left(ke\right)cosk\tau \end{array}\begin{array}{l}{S}_n^{\prime }=\sum _{k=1}^{\infty }{k}^n{J}_k^{\prime }\left(ke\right)sink\tau \\ C_n^{\prime }=\sum _{k=1}^{\infty }{k}^n{J}_k^{\prime }\left(ke\right)cosk\tau \end{array}$

where a prime indicates the series is constructed from derivatives of Bessel functions and $\tau =\omega t$ as before. Derivatives with respect to this mean anomaly increase the index while changing the circular function:

$\begin{array}{l}\frac{\partial S_n}{\partial \tau }=C_{n+1}\\ \frac{\partial C_n}{\partial \tau }=-S_{n+1}\end{array}\begin{array}{l}\frac{\partial S_n^{\prime }}{\partial \tau }=C_{n+1}^{\prime }\\ \frac{\partial C_n^{\prime }}{\partial \tau }=-S_{n+1}^{\prime }\end{array}$

Derivatives with respect to eccentricity increase the index and alternate whether the series is constructed from derivatives or not:

$\begin{array}{l}\frac{\partial S_n}{\partial e}=S_{n+1}^{\prime }\\ \frac{\partial C_n}{\partial \tau }=C_{n+1}^{\prime }\end{array}\begin{array}{l}\frac{\partial \left(e{S}_n^{\prime }\right)}{\partial e}=\frac{1-e^2}{e}{S}_{n+1}\\ \frac{\partial \left(e{C}_n^{\prime }\right)}{\partial e}=\frac{1-e^2}{e}{C}_{n+1}\end{array}$

The defining differential equation for Bessel functions has been used here to replace second derivatives.

In terms of this notation, simple dynamical variables are

$\begin{array}{lll}cosu=-\frac{e}{2}+2C_{-1}^{\prime }& & sinu=\frac{2}{e}{S}_{-1}\\ \\ \frac{r}{a}=1+\frac{e^2}{2}-2e{C}_{-1}^{\prime }& & \frac{a}{r}=1+2C_0\\ \\ \frac{x}{a}=-\frac{3e}{2}+2C_{-1}^{\prime }& & \frac{y}{a}=\frac{2\sqrt{1-e^2}}{e}{S}_{-1}\\ \\ cos\phi =-e+\frac{2(1-e^2)}{e}{C}_0& & sin\phi =2\sqrt{1-e^2}{S}_0^{\prime }\end{array}$

and their combinations are

$\begin{array}{l}{cos}^{2}u=\frac{1}{2}-\frac{4}{e^2}{C}_{-2}+\frac{4}{e}{C}_{-1}^{\prime }\\ \\ {sin}^{2}u=\frac{1}{2}+\frac{4}{e^2}{C}_{-2}-\frac{4}{e}{C}_{-1}^{\prime }\\ \\ sinucos\phi =-\frac{2}{e^2}{S}_{-1}+\frac{2(1-e^2)}{e}{S}_0^{\prime }\\ \\ sinusin\phi =\sqrt{1-e^2}[\frac{1}{2}-\frac{2(1-e^2)}{e^2}{C}_0+\frac{2}{e}{C}_{-1}^{\prime }]\\ \\ \frac{a^2}{r^2}cos\phi =2C_1^{\prime }\\ \\ \frac{a^2}{r^2}sin\phi =\frac{2\sqrt{1-e^2}}{e}{S}_1\end{array}$

Algebraic relationships among the series can be determined from the algebraic relationships of these variables. For example, multiplying the radial variable and its inverse gives

$\begin{array}{c}\frac{r}{a}·\frac{a}{r}=(1+\frac{e^2}{2}-2e{C}_{-1}^{\prime })(1+2C_0)\equiv 1\\ \\ C_{-1}^{\prime }{C}_0=\frac{e}{8}-\frac{1}{2}{C}_{-1}^{\prime }+\frac{2+e^2}{4e}{C}_0\end{array}$

Equating the expressions for $\frac{r}{a}cos\phi$ and $\frac{x}{a}$ produces the same result. Using expressions for the y-coordinate one has

$\begin{array}{c}\frac{r}{a}sin\phi =(1+\frac{e^2}{2}-2e{C}_{-1}^{\prime })·2\sqrt{1-e^2}{S}_0^{\prime }=\frac{2\sqrt{1-e^2}}{e}{S}_{-1}=\frac{y}{a}\\ \\ C_{-1}^{\prime }{S}_0^{\prime }=-\frac{1}{2e}{S}_{-1}+\frac{2+e^2}{4e}{S}_0^{\prime }\end{array}$

This relationship arises from adding the squares of circular functions,

$\begin{array}{c}{sin}^2\phi +{cos}^2\phi =4(1-e^2){S_0^{\prime }}^2+(-e+\frac{2(1-e^2)}{e}{C}_0{)}^2\equiv 1\\ \\ (1-e^2)C_0^2+(e{S}_0^{\prime }{)}^2=\frac{e^2}{4}+e^2{C}_0\end{array}$

but there does not appear to be a simple way to determine each of these squared series separately. This can be seen by beginning to evaluate the expansion for the squared inverse of the radial variable, which would be needed for a separate expression for $C_0^2$ . The constant term is a known integral,

$\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}\frac{a^2}{r^2}d\tau =\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}\frac{du}{1-ecosu}=\frac{1}{\sqrt{1-e^2}}$

but the remaining coefficients are not so simple:

$\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}\frac{a^2}{r^2}cosk\tau d\tau =\frac{2}{\pi }\underset{0}{\overset{\pi }{\int }}\frac{cos\left[k\right(u-esinu\left)\right]}{1-ecosu}du$

The denominator can be expanded as an infinite geometric series and rearranged to terms of the form $cosmu$ for explicit evaluation, but that means that each of these coefficients is itself an infinite series. In fact this same infinite series appears for the coefficients of the expansion of φ mentioned above but not evaluated.

There are clearly some simplifications among products of these interesting series that are akin to the product formulae of circular functions. To be continued...

Uploaded 2018.08.30 — Updated 2020.11.15 analyticphysics.com